# RBC的欧拉方程迭代
rm(list = ls())
library(ggplot2)
library(magrittr)
library(signal)
devtools::load_all()

# 参数设置
a <- seq(0.06,6.5,0.06)
beta <- .98
delta <-  .1
theta <- .36
numits <- length(a)

#----------Carroll (2006)的内生格点法-------
# 函数定义: u = log(cc), 对数效用下偏导和偏导的逆是一样的
NP <- UP <- function(cc) 1/cc

# 初始化
m <- cc <- seq(1,length(a)) %>% matrix(ncol = 1)

for (i in 1:numits) {
  # 插值计算c(t+1),10式中mt的定义
  m_next <- a^theta + (1-delta)*a
  c_next <- interp1(m[,ncol(m)],cc[,ncol(cc)],m_next,extrap = T)

  # 笔记中的13式反算出ct
  cc <- NP(beta * (theta * a^(theta-1)+1-delta) * UP(c_next)) %>% cbind(cc,.)

  # 更新m
  m <- (a + cc[,ncol(cc)]) %>% cbind(m,.)
}

# 通过m反算出k
mfun <- function(kt,mt) kt^theta+(1-delta)*kt-mt
k <- numeric(nrow(m))
for (i in 1:numits) {
  k[i] <- uniroot(mfun,c(0,60),mt = m[i,ncol(m)])$root
}

# 画图
ggplot(data.frame(kt = k, knext = a),aes(x = kt,y = knext)) + geom_line()
data.frame(k = k, cc = cc[,ncol(cc)], m = m[,ncol(m)],knext = a)

# 模拟时间序列
# 因为欧拉方程需要t+1期的c和t+1期的k才能算出t期的c,所以只能反向递归，从终值往前
numT <- 30
kT <- cT <- numeric(numT)
cT[length(cT)] <- 1.25
for (i in numT:2) {
  kT[i] <- knext <- interp1(cc[,ncol(cc)], a, cT[i], extrap = T)
  cT[i-1] <- NP(beta * UP(cT[i]) * (theta * knext^(theta-1)+1-delta))
}
cT
kT

# ---------标准解法--------
# 初始化
k <- cc <- seq(1,10,0.06) %>% matrix(ncol = 1)

for (i in 1:numits) {
  knext <- k^theta + (1-delta)*k-cc[,ncol(cc)]
  cnext <- interp1(k,cc[,ncol(cc)],knext, extrap = T)
  cc <- (1/(beta/cnext * (theta * knext^(theta - 1)+1-delta))) %>% cbind(cc,.)
}
picdata <- data.frame(k = k, cc = cc[,ncol(cc)])
picdata$knext <- picdata$k^theta+(1-delta)*k-picdata$cc
ggplot(picdata,aes(x = k,y = knext)) + geom_line()


#---------相图---------------
k <- seq(0,20,0.06)
cc <- k^theta-delta*k
phase <- data.frame(k = k, cc = cc)
kvt <- ((1/beta-1+delta)/theta)^(1/(theta-1))

ggplot(phase, aes(x = k, y = cc)) + geom_line(size = 1.5, color = I('steelblue1')) +
  geom_vline(xintercept = kvt,size = 1.5, color = I('steelblue1')) + geom_line(data = picdata, aes(x = k, y = cc))+
  annotate('segment', x = picdata$k[c(30,45)], xend = picdata$k[c(40,55)],
           y = picdata$cc[c(30,45)],yend = picdata$cc[c(40,55)],
           arrow = arrow(length = unit(.2,'cm')))+
  annotate('segment', x = picdata$k[c(130,115)], xend = picdata$k[c(120,105)],
           y = picdata$cc[c(130,115)],yend = picdata$cc[c(120,105)],
           arrow = arrow(length = unit(.2,'cm'))) + labs(y = 'c')+
  annotate('segment', x = 4, xend = 9,
           y = 0.5,yend = 0.5,
           arrow = arrow(length = unit(.2,'cm'))) + labs(y = 'c')+
  annotate('segment', x = 9, xend = 4,
           y = 1.6,yend = 1.6,
           arrow = arrow(length = unit(.2,'cm'))) + labs(y = 'c')+
  annotate('segment', x = 10, xend = 10,
           y = 1.5,yend = 1,
           arrow = arrow(length = unit(.2,'cm'))) + labs(y = 'c')+
  annotate('segment', x = 2.5, xend = 2.5,
           y = 1,yend = 1.5,
           arrow = arrow(length = unit(.2,'cm'))) + labs(y = 'c')+
  theme_bw()
# ggsave('data-raw/phase.pdf')

#---------使用dynare对数线性化求解------
library(DynareR)
library(matlabr)
library(R.matlab)
mtl_src <- '
cd ./data-raw
dynare rbc
'
run_matlab_code(mtl_src)
ans <- readMat('data-raw/rbc_results.mat')
ans$oo.[[4]]

#---------使用deSolve包--------
library(deSolve)
library(bvpSolve)
nroot <- function(x,n) abs(x)^n*sign(x)

rbc <- function(tm,y,parms){
  with(as.list(y),{
    dk <- nroot(K,theta)-delta*K-ct
    dc <- beta*ct*(theta*nroot(nroot(K,theta)-delta*K-ct+K,theta-1)+1-delta) -ct
    list(c(dk,dc))
  })
}
ode(y = c(K=5.54,ct=1.23),times = seq(0,5,0.2),func = rbc)
bvpshoot(yini = c(K=5.54,ct = NA),yend = c(K = NA,ct = 1.3),func = rbc, x = seq(0,5,0.2))
